set.seed(5168)
## responsibility attribution ##

	library(foreign)
	library(mnormt)
	library(ordinal)

  data <- read.dta('Denmark (2014).dta')

	model <- clmm(as.factor(Resp_attribution) ~ DK_V2_Prime_minister + DK_V2_Cabinet + DK_V2_Support + DK_V2_Opposition + DK_Perc_seats + DK_V2_Prime_minister * DK_Perc_seats + DK_V2_Cabinet * DK_Perc_seats + DK_V2_Support * DK_Perc_seats + DK_V2_Opposition*DK_Perc_seats + (1 | Party) + (1 | pid), data = data, HESS = TRUE, link = 'probit')
	summary(model)

	pars <- as.matrix(model$coefficients)
	vce <- vcov(model)
	vce <- vce[-c(14:15), ]      
	vce <- vce[, -c(14:15)]      
	
	simCoef <- rmnorm(1000, as.vector(pars),as.matrix(vce))

	summary(model)
	summary(simCoef)

	pmHi <- c()
	pmLo <- c()
	pmMed <- c()

	partHi <- c()
	partLo <- c()
	partMed <- c()

	supHi <- c()
	supLo <- c()
	supMed <- c()

	oppHi <- c()
	oppLo <- c()
	oppMed <- c()

  otherHi <- c()
  otherLo <- c()
  otherMed <- c()


	for(i in 0:50){
	## pm party 

		agg <- simCoef[, 5] * 1 + 
				simCoef[, 6] * 0 + 
				simCoef[, 7] * 0 + 
				simCoef[, 8] * 0 + 
				simCoef[, 9] * i + 
				simCoef[, 10] * i + 
				simCoef[, 11] * 0 + 
				simCoef[, 12] * 0 + 
				simCoef[, 13] * 0 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		pmHi[i + 1] <- quantile(pred, 0.975)
		pmLo[i + 1] <- quantile(pred, 0.025)
		pmMed[i + 1] <- quantile(pred, 0.5)


	## partner party
	agg <- simCoef[, 5] * 0 + 
	  simCoef[, 6] * 1 + 
	  simCoef[, 7] * 0 + 
	  simCoef[, 8] * 0 + 
	  simCoef[, 9] * i + 
	  simCoef[, 10] * 0 + 
	  simCoef[, 11] * i + 
	  simCoef[, 12] * 0 + 
	  simCoef[, 13] * 0 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		partHi[i + 1] <- quantile(pred, 0.975)
		partLo[i + 1] <- quantile(pred, 0.025)
		partMed[i + 1] <- quantile(pred, 0.5)

	## support party
	agg <- simCoef[, 5] * 0 + 
	  simCoef[, 6] * 0 + 
	  simCoef[, 7] * 1 + 
	  simCoef[, 8] * 0 + 
	  simCoef[, 9] * i + 
	  simCoef[, 10] * 0 + 
	  simCoef[, 11] * 0 + 
	  simCoef[, 12] * i + 
	  simCoef[, 13] * 0 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		supHi[i + 1] <- quantile(pred, 0.975)
	  supLo[i + 1] <- quantile(pred, 0.025)
	  supMed[i + 1] <- quantile(pred, 0.5)

	## opp party 
  	agg <- simCoef[, 5] * 0 + 
  	  simCoef[, 6] * 0 + 
  	  simCoef[, 7] * 0 + 
  	  simCoef[, 8] * 1 + 
  	  simCoef[, 9] * i + 
  	  simCoef[, 10] * 0 + 
  	  simCoef[, 11] * 0 + 
  	  simCoef[, 12] * 0 + 
  	  simCoef[, 13] * i 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		oppHi[i + 1] <- quantile(pred, 0.975)
		oppLo[i + 1] <- quantile(pred, 0.025)
		oppMed[i + 1] <- quantile(pred, 0.5)

    ## No seats
    agg <- simCoef[, 5] * 0 + 
      simCoef[, 6] * 0 + 
      simCoef[, 7] * 0 + 
      simCoef[, 8] * 0 + 
      simCoef[, 9] * i + 
      simCoef[, 10] * 0 + 
      simCoef[, 11] * 0 + 
      simCoef[, 12] * 0 + 
      simCoef[, 13] * 0 
    
    pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
    pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
    pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
    pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
    pred5 <- 1 - pred1 - pred2 - pred3 - pred4
    
    pred <- pred5 * 5 + 
      pred4 * 4 + 
      pred3 * 3 + 
      pred2 * 2 + 
      pred1
    
    otherHi[i + 1] <- quantile(pred, 0.975)
    otherLo[i + 1] <- quantile(pred, 0.025)
    otherMed[i + 1] <- quantile(pred, 0.5)
    
    }


## now paint a pretty picture.... ##

	seats <- seq(0, 50, 1)       
  	
	#Now create the figure for Denmark
	par(mar=c(3.6,4,1,2))    
	
	plot(1, 1, type = 'n',
		xlim = c(0, 50),
		ylim = c(1, 5),
		xlab = '',
		ylab = '',
		main = '',
		axes = FALSE)
	axis(1,at=seq(0,50,25),labels=T)   
	axis(2)

lines(seats, pmMed, col = rgb(1, 0, 0, 0.2), lwd = 3)

## darken bands in high density range ##
	seatHi <- quantile(data$DK_Perc_seats[data$DK_V2_Prime_minister == 1], 0.9, na.rm = TRUE)
	seatLo <- quantile(data$DK_Perc_seats[data$DK_V2_Prime_minister == 1], 0.1, na.rm = TRUE)
	lines(seats[seats >= seatLo & seats <= seatHi], pmMed[seats >= seatLo & seats <= seatHi], col = rgb(1, 0, 0, 0.3), lwd = 15)
	
	lines(seats, partMed, col = rgb(0, 0, 1, 0.2), lwd = 3)

	seatHi <- quantile(data$DK_Perc_seats[data$DK_V2_Cabinet == 1], 0.9, na.rm = TRUE)
	seatLo <- quantile(data$DK_Perc_seats[data$DK_V2_Cabinet == 1], 0.1, na.rm = TRUE)
	lines(seats[seats >= seatLo & seats <= seatHi], partMed[seats >= seatLo & seats <= seatHi], col = rgb(0, 0, 1, 0.3), lwd = 15)

	lines(seats, supMed, col = rgb(0, 1, 0, 0.2), lwd = 3)

	seatHi <- quantile(data$DK_Perc_seats[data$DK_V2_Support == 1], 0.9, na.rm = TRUE)
	seatLo <- quantile(data$DK_Perc_seats[data$DK_V2_Support == 1], 0.1, na.rm = TRUE)
	lines(seats[seats >= seatLo & seats <= seatHi], supMed[seats >= seatLo & seats <= seatHi], col = rgb(0, 1, 0, 0.3), lwd = 15)

	lines(seats, oppMed, col = rgb(0, 0, 0, 0.2), lwd = 3)

	seatHi <- quantile(data$DK_Perc_seats[data$DK_V2_Opposition == 1], 0.9, na.rm = TRUE)
	seatLo <- quantile(data$DK_Perc_seats[data$DK_V2_Opposition == 1], 0.1, na.rm = TRUE)
	lines(seats[seats >= seatLo & seats <= seatHi], oppMed[seats >= seatLo & seats <= seatHi], col = rgb(0, 0, 0, 0.3), lwd = 15)

  points(0, y = otherMed[1], col = rgb(1, 0, 1, 0.3), lwd = 15) 

	mtext("Responsibility attribution", side=2, line=2.5, cex=1.5)
	mtext("Perceived seat %", side=1, line=2.2, cex=1.5)           


############
###LEGEND###
############

legend('bottomright', legend = c('PM', 'Partner', 'Support', 'Opposition', 'No seats'), fill = c(rgb(1, 0, 0, 0.3),	
	rgb(0, 0, 1, 0.3), rgb(0, 1, 0, 0.3), rgb(0, 0, 0, 0.3), rgb(1, 0, 1, 0.3)), bty = 'n', border = NA, , cex=1.25)

###################
###EXPORT FIGURE###
###################

dev.copy(png,'Figure 5b DK 2014 - party types over seat share.png', height = 490, width = 633)	
dev.off()  
